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SOLUTION  TO  THE  COMPRESSIBLE  NAVIER-STOKES  EQUATIONS 
OF  MOTION  BY  CHEBYSHEV  POLYNOMIALS  FOR 
LAMINAR  SHOCK-BOUNDARY  LAYER  FLOW 


1.  INTRODUCTION 

Gottlieb  et.  al.  (reference  1)  applied  pseudospectral  methods  to  the  solution  of  the 
one  dimensional  propagating  shock  wave  problem.  They  utilized  a  low  pass  spectral 
filter  which  they  developed  together  with  a  Shuman  filter,  applied  to  the  flow  on  either 
side  of  the  shock  wave  but  not  across  the  shock  front  itself.  The  shock  location  was 
determined  by  examination  of  the  spectral  coefficients.  Since  then,  the  present  author 
has  developed  new  techniques  for  use  with  pseudospectral  methods  which  have  greatly 
increased  their  utility  in  solving  inviscid  flows  with  single  or  multiple  discontinuities. 
Pseudospectral  methods  have  been  used  by  the  author  to  solve  many  classes  of  com¬ 
plicated  time  dependent  compressible  flows  using  the  full  Euler  equations  of  motion 
(references  2  through  6).  Results  have  shown  that  flow  discontinuities  such  as  shock 
waves  or  contact  surfaces  are  properly  resolved  as  sharp  discontinuities.  Solutions  for 
transonic  airfoil  flows  at  subcritical  and ‘supercritical  conditions  (reference  6)  were  ob¬ 
tained  more  recently  and  proved  that  full  pseudospectral  computational  methods  could 
also  successfully  treat  compressible  inviscid  flows  about  non-planar  geometries.  With 
the  completion  of  the  airfoil  work,  the  author  has  shown  that  pseodospectral  compu¬ 
tational  methods  are  fully  suitable  for  solving  any  class  of  inviscid.  time  dependent, 
compressible  flow  using  the  Euler  equations  of  motion. 

The  next  logical  step  is  to  turn  to  the  full  viscous  equations  of  motion  namely,  the 
Navier-Stokes  equations.  Orszag  is  the  most  notable  in  the  field  when  one  considers 
stability  and  transition  analyses  of  incompressible  hydrodynamic  boundary  layer  and 
channel  flows  (reference  7  for  example).  He  was  the  first  to  apply  pseudospectral 
methods,  treating  flow  stability  problems  for  laminar  hydrodynamic  boundary  layers. 
However,  no  work  has  as  yet  been  done  for  compressible  flows.  Even  when  one  considers 
the  solution  to  compressible  external  flows,  pseudospectral  methods  have  not  been  used 
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at  ail  due  to  difficulties  in  resolving  discontinuities.  The  present  work  remedies  this 
difficulty  and  shows  that  when  techniques  developed  by  the  author  are  implemented,  full 
pseudospectral  computational  methods  can  successfully  solve  viscous,  time  dependent 
compressible  flows  with  multiple  discontinuities  and  flow  separation. 


2.  GOVERNING  EQUATIONS 

The  full  two-dimensional,  time  dependent,  compressible  Navier-Stokes  equations 
of  motion  cast  in  conservation  law  form  are  shown  below. 
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where  U  ,  E  and  F  are  vectors  whose  elements  are: 
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where  crx  ,cry  are  the  normal  stresses, rxy  and  ryx  are  the  shear  stresses,  A  and  \i  are 
viscosity  coefficients  (Sutherland’s  relation  is  used  since  the  present  work  deals  with 
laminar  flow)  and  k  is  the  coefficient  of  thermal  conductivity.  The  pressure  is  obtained 
from  the  following 

e  =——-  +  ]■  p{u2 +  v2),  (2e) 

(7-1)  2 

The  physical  flow  variables  are  non-dimensionalized  in  the  following  manner. 
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The  normal  and  shear  stress  terms  are  non-dimensionalized  by  the  free  stream 
pressure  head,  (piUf)  ■  Subscript  one  denotes  free  stream  properties  upstream  of  the 
incident  shock  wave.  With  respect  to  non-dimensional  flow  variables,  equation  one 
becomes: 
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with  the  Prandtl  number  Pr  defined  bv, 


Pr  = 


Tins  completes  the  non-dimensionalization  of  the  physical  flow  variables  and  the 
conversion  of  the  Navier  Stokes  equations  to  non-dimensional  physical  flow  variables. 
However,  it  still  remains  to  transform  these  equations  into  a  suitable  computational 
spare.  This  is  discussed  below.  Several  coordinate  transformations  are  applied  to 
generate  an  appropriate  distribution  of  points  in  the  flow  field.  Appropriate  here  means 
many  points  in  regions  of  large  gradients  and  simultaneously  few  points  in  regions  of 
small  gradients.  The  final  computational  coordinates  are  obtained  using  a  sequence  of 
four  coordinate  transformations.  Namelv. 
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The  terms  Ci,  C3,  ,4x,,4y  and  ,  a  are  transformation  clustering  constants  which 
affect  the  distribution  of  points  in  the  computational  and  physical  domains.  The  final 
form  of  the  Navier  Stokes  equations  becomes. 
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All  spatial  derivatives  appearing  in  equation  8  are  calculated  by  pseudospectral 
means.  In  the  present  work  this  involves  the  use  of  chebvshev  polynomials.  The  time 
derivative  Ut  appearing  in  equation  8  is  evaluated  using  finite  differences.  Specifically, 
the  Adams  Bashforth  algorithm  is  used.  The  resultant  difference  form  of  equation  S  is 
given  by. 
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The  term  D,,;  is  an  artificial  viscosity.  In  the  present  work  a  fourth  order  artifi¬ 
cial  viscosity  is  utilized  and  is  computed  using  finite  differences.  The  finite  difference 
representation  is  shown  below. 
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The  terms  Dz  and  Dy  are  smoothing  constants. 


3.  PSEUDOSPECTRAL  METHODS 

Pseudospectral  solution  techniques  involve  the  use  of  series  of  functions  to  represent 
the  global  properties  of  a  flow  field  and  its  spatial  derivatives.  In  the  present  work 
Chebyshev  polynomials  are  used.  They  are  represented  by  Tn(x)  where 


Tn(x)  =  cos  [narccos(x)] 
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or 
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A  function  of  a  single  spatial  variable  and  time  such  as  F(x.t)  may  be  represented 
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n  =  0 

The  time  dependence  is  represented  entirely  in  the  spectral  coefficients  A„(t)  while, 
the  spatial  dependence  is  represented  in  the  Chebyshev  polynomials  T„(.r).  The  Cheby¬ 
shev  polynomials  are  evaluated  at  discrete  points  xj  where 


,r 


j 


cos 


(14) 


where  Nz  is  the  total  number  of  modes  used  to  represent  the  spatial  variation  of 
the  function  F(x,t).  The  spatial  derivative  of  the  function  F(x,t)  is  represented  as 
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The  An's  are  determined  from  equation  thirteen.  Inverse  FFT’s  are  used  to  obtain 
the  An' s  from  the  known  functional  values  F(x,t)  at  the  known  collocation  points  Xj. 
The  spectral  coefficients  of  the  spatial  derivative  ,  ,  are  determined  from  the 

recurrence  relation  equation  15b.  Direct  FFT's  are  used  to  evaluate  the  sum  in  equation 
thirteen  to  obtain  the  functional  values  at  t+<5t.  The  low  pass  spectral  filter  developed 
by  Gottlieb  at  ICASE  is  used  to  damp  spectral  oscillations.  It  is  shown  below. 
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where  K  is  the  spectral  wavenumber  and  Kmax  is  the  maximum  wavenumber 
corresponding  to  the  total  number  of  collocation  points. 
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4.  RESULTS 

The  following  free  stream  conditions  were  selected  to  match  one  of  the  wind  tunnel 
experiments  reported  in  reference  S.  An  experimentally  measured  surface  pressure 
distribution  is  included  in  the  above  reference  so  that  a  direct  comparison  can  be  made 
between  the  experimental  and  numerical  distributions.  The  free  stream  mach  and 
unit  Reynolds  numbers  are  2.05  and  095,000  per  foot  respectively.  The  free  stream 
Prandtl  number  is  0.71.  The  incident  shock  wave  was  generated  by  a  six  degree  wedge. 
The  physical  extent  of  the  computational  domain  is  0  <  ;/  <  0.3  foot  and  —0.2  < 
x  <  +0.2.  No  attempt  was  made  to  select  an  optimum  shape  for  the  computational 
boundary  ,  such  as  alligning  constant  coordinate  lines  with  the  incident  and  reflected 
shock  wave  system,  since  the  purpose  of  the  present  work  is  to  determine  whether 
pseudospectral  computational  techniques  can  actually  solve  the  full  time  dependent, 
compressible  Navier-Stokes  equations  of  motion  with  discontinuities.  The  flat  plate 
surface  lies  in  the  range  —.14  <  x  <  +0.2  so  that  several  grid  points  lie  ahead  of  the 
plate  leading  edge.  Sixty  four  modes  each  were  used  to  represent  the  flow  field  in  the 
x  and  y  directions.  The  physical  space  computational  boundary  is  shown  in  Figure 
1  along  with  constant  coordinate  lines  showing  the  degree  of  clustering  in  both  the 
x  and  y  directions.  As  can  be  seen,  points  are  highly  clustered  in  the  y  direction  at 
the  surface.  This  was  done  to  properly  resolve  the  flow  separation  zone  and  attendent 
shock  structure  which  arises  very  close  to  the  plate  surface.  Points  are  only  mildly 
clustered  along  the  x-direction  since  an  airfoil  coordinate  transformation  which  clusters 
points  about  the  leading  and  trailing  edges  was  used  with  the  clustering  parameters 
detuned.  The  inflow  boundary  conditions  were  to  keep  all  flow  variables  fixed  because 
of  supersonic  inflow.  Along  the  upper  boundary  flow  variables  wore  held  fixed  at 
post  incident  shock  conditions  since  for  the  flow  considered  the  reflected  shook  does 
not  extend  to  the  upper  boundary  location.  At  the  outflow  boundary  conditions  at 
each  point  were  set  to  those  at  the  next  upstream  point  (zero  th  order  extrapolation) 
throughout  the  calculation.  Along  the  bottom  of  the  computational  boundary  either 
a  plane  of  symmetry  or  wall  surface  was  present.  Reflective  boundary  conditions  were 
used  at  points  ahead  plate  leading  edge  namely.  t/,,i  =  a, ,-j.  On  the  plate  surface. u,A  = 
—  u,,?.  Also,  along  the  entire  bottom  boundary  e,.i  =-  e,  _■».  p,.\  —  P>.2  and  e,  \  =  e,:i. 
with  (  denoting  specific  internal  energy.  This  last  condition  being  applied  along  the 
plate  surface  since  the  case  being  treated  is  <m  adiabatic  wall.  The  time  step  size  was 
determined  from  the  following. 
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For  results  presented  herein  the  courant  number  CN  was  held  fixed  at  2.5.  No 
attempt  was  made  to  replace  A£mtn  or  A(min  with  the  actual  respective  values  at  the 
points  of  maximum  eigenvalue.  The  flowfield  was  initialized  in  two  steps.  First,  the 
exact  euler  shock  structure  (  incident  and  reflected  )  including  post  shock  zones  was  put 
in.  Figure  2  shows  the  initial  shock  system.  This  system  is  overlayed  on  the  constant 
grid  lines  plot  to  show  the  relative  grid  resolution  over  the  entire  shock  system  in  Figure 
3.  The  second  step  was  to  put  in  a  boundary  layer  field  using  a  blasius  velocity  profile 
and  the  Crocco  relation  to  obtain  the  density  profile  at  constant,  pressure  at  each  point 
along  the  plate  surface.  The  pseudospectral  code  used  to  obtain  the  present  results  was 
run  on  the  NRL  CRAY  XMP/12.  The  code  takes  0.75  cpu  second  per  iteration  at  a  grid 
resolution  of  64  by  64.  The  solution  converged  in  25,000  iterations  or  a  little  over  five 
cpu  hours.  No  optimization  of  any  kind  was  employed  since  the  purpose  of  this  work 
was  to  determine  if  pseudospectral  techniques  could  successfully  solve  the  compressible 
Navier-Stokes  equations.  Pressure  contours  of  the  entire  computational  field  are  shown 
in  Figure  4.  There  are  no  oscillations  in  the  field.  The  incident  shock  is  slightly  smeared 
for  several  reasons.  First,  the  grid  system  is  not  shock  alligned.  Second,  the  x-direction 
grid  resolution  is  only  moderate.  Finally,  at  a  courant  number  of  2.5  the  magnitude 
of  the  artificial  viscosity  coefficients  required  for  stability  is  relatively  large  to  that 
required  to  keep  the  shock  waves  stable.  This  last  condition  is  primarily  responsible  fro 
the  shock  smearing.  With  a  courant  number  of  0.7  the  thickness  of  the  shock  is  only 
one  third  of  that  shown  in  the  figures  herein  included.  However,  the  penalty  that  one 
incurs  for  the  reduction  of  shock  smearing  is  an  increase  on  CPU  time.  For  the  courant 
numbers  thus  mentioned,  the  CPU  time  would  be  increased  by  almost  a  factor  of  four. 
It  was  deemed  more  important  for  the  present  work  to  obtain  a  converged  solution 
while  minimizing  the  CPU  time.  An  alternative  would  be  to  run  at  maximum  courant 
number  untill  convergence  is  almost  reached  and  then  reduce  the  courant  number  and 
dissipation  constants  to  have  the  best  of  both  worlds.  The  reflected  shock  wave  has 
split  into  two  compression  fans.  One  is  well  above  the  separation  zone,  while  the  second 
is  the  re-compression  fan  at  the  point  of  re-attachment  of  the  separated  shear  layer  on 
the  plate  surface.  A  blowup  of  the  wall  region  pressure  contours  (0  <  y  <  0.0S fool)  is 
shown  in  Figure  5.  Corresponding  velocity  vector  and  mass  flow  per  unit  area  contours 
are  shown  in  Figures  6  and  7.  The  separation  zone  is  clearly  evident.  Comparison 
of  the  numerical  and  experimental  surface  pressure  distributions  can  be  made  from 
Figure  S.  The  starting  location  of  the  separation  zone  is  in  good  agreement  with 


the  data.  The  experimentally  measured  plateau  pressure  non-dimensionalized  by  pre¬ 
impingement  free  stream  pressure  is  1.25.  The  numerical  value  is  about  1.29  or  a  little 
more  than  three  per  cent  higher.  The  re-compression  point  location  is  in  very  good 
agreement.  The  numerical  solution  value  is  0.155  foot  from  the  plate  leading  edge 
while  the  experimental  value  is  0.1G5  foot.  The  length  of  the  separation  zone  from 
the  present  work  is  0.09  foot  or  1.0S  inch.  The  actual  length  is  about  one  inch  as 
determined  from  a  schlieren  photograph.  The  height  of  the  separation  zone  is  about 
the  same  as  the  experimental  value.  (This  is  all  that  can  be  said  since  the  schlieren 
photograph  in  reference  S  is  too  small  to  measure  the  separation  height. )  The  numerical 
pressure  distribution  through  the  recompression  point  is  more  rounded  than  the  sharp 
discontinuity  of  the  experimentally  obtained  distribution.  The  sharpness  is  in  part  due 
to  the  larger  spacing  of  the  experimental  data  points  compared  to  the  numerical  grid 
point  spacing.  Thereafter,  the  numerical  solution  is  in  excellent  agreement  with  the 
experimental  data,  falling  right  on  top  of  the  experimental  data  points. 

Vertical  profile  plots  of  density,  u-velocitv.  energy  and,  mass  How  per  unit  area 
are  shown  in  Figures  9  through  12  for  the  full  computational  domain  as  well  as  for  the 
wall  region.  The  vertical  extent  is  0.30  and  0.01  foot  respectively,  while  the  horizontal 
range  of  the  plots  covers  each  of  the  x  locations.  The  horizontal  spacing  of  the  lines  is 
proportional  to  the  grid  spacing  in  the  x-dircction  and  is  not  equal  to  it.  The  change  in 
horizontal  location  of  any  curve  from  its  surface  position  is  representative  of  the  change 
in  the  magnitude  of  the  flow  variable  being  plotted  from  the  surface  value.  There  are 
no  numerical  oscillations  present  at  the  wall.  The  separation  zone  is  dearly  shown  in 
Figure  10. 


5.  CONCLUSIONS 

The  present  work  has  shown  that  pseudospectnd  computational  techniques  are 
indeed  able  to  accurately  solve  the  compressible  Navier-Stokes  equations  for  flows  si¬ 
multaneously  including  separation  zones  and  shock  waves.  No  numerical  oscillations 
are  present  in  the  solution.  It  still  remains  to  incorporate  existing  acceleration  tech¬ 
niques  as  well  as  to  develop  new  ones  to  reduce  the  large  computer  time  requirt'd  to 
obtain  a  converged  solution  to  a  more  practical  level.  Work  is  currently  proceeding 
along  these  lines. 
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